Community differentiation of rhizosphere microorganisms and their responses to environmental factors at different development stages of medicinal plant Glehnia littoralis

Rhizosphere microorganisms play a key role in affecting plant quality and productivity through its interaction with plant root system. To figure out the bottleneck of the decline of yield and quality in the traditional Chinese medicinal herbs Glehnia littoralis they now encounter, it is important to study the dynamics of rhizosphere microbiota during the cultivation of G. littoralis. In the present study, the composition, diversity and function of rhizosphere microbes at different development stages of G. littoralis, as well as the correlation between rhizosphere microbes and environmental factors were systematically studied by high-throughput sequencing. There were significant differences between the rhizosphere microbes at early and middle-late development stages. More beneficial bacteria, such as Proteobacteria, and more symbiotic and saprophytic fungi were observed at the middle-late development stage of G. littoralis, while beneficial bacteria such as Actinobacteria and polytrophic transitional fungi were abundant at all development stages. The results of redundancy analysis show that eight environmental factors drive the changes of microflora at different development stages. pH, soil organic matter (SOM) and available phosphorus (AP) had important positive effects on the bacterial and fungal communities at the early development stage; saccharase (SC) and nitrate nitrogen (NN) showed significant positive effects on the bacterial and fungal communities at the middle and late stages; while urease (UE), available potassium (AK), and alkaline phosphatase (AKP) have different effects on bacterial and fungal communities at different development stages. Random forest analysis identified 47 bacterial markers and 22 fungal markers that could be used to distinguish G. littoralis at different development stages. Network analysis showed that the rhizosphere microbes formed a complex mutualistic symbiosis network, which is beneficial to the growth and development of G. littoralis. These results suggest that host development stage and environmental factors have profound influence on the composition, diversity, community structure and function of plant rhizosphere microorganisms. This study provides a reference for optimizing the cultivation of G. littoralis.


INTRODUCTION
The plant root system, rhizosphere microorganisms and rhizosphere soil constitute a plant rhizosphere microecosystem. In this microecosystem, biotic factors (plant genotypes, plant developmental stages, invasive pathogenic microorganisms) and abiotic factors (soil composition, soil management and climatic conditions) in soil influence the composition, diversity, structure and function of rhizosphere microbial communities (Vacheron et al., 2013;Gouda et al., 2018;Li et al., 2020;Ma et al., 2021). In turn, rhizosphere microorganisms influence plant root exudates and even plant metabolome (Badri et al., 2013), which involves improving plant physiology and resistance to pathogens through various mechanisms (Zakry et al., 2012;Bai et al., 2022), affecting plant health and growth (Trivedi et al., 2020;Xiong et al., 2021). Rhizosphere microorganisms also affect soil evolution and play a vital role in the conversion of poor soil and low-quality soil into cultivable soil (Gouda et al., 2018). In summary, multitrophic interactions between plants, microorganisms and environmental factors lead to the formation of complex symbiotic networks in rhizosphere microecosystem. This symbiotic network dynamically affects rhizosphere communities and alters plant phenotypes (Lu et al., 2018).
Beneficial rhizosphere microorganisms can enhance plant roots' vitality, promote plant growth, increase plant yield and improve resistance to phytopathogens (Razavi et al., 2016). Conversely, harmful microorganisms can lead to plant disease, the inhibition of root growth, the suppression of plant growth, and crop failure (Kremer, 2006). The study of the symbiosis and interactions between medicinal plants, environmental factors and microbes can help regulate the microbiota of medicinal plants themselves and their surroundings, contribute to the goal of high quality and high yield of medicinal plants (Karthikeyan et al., 2008;Khamna, Yokota & Lumyong, 2009;Köberl et al., 2013;Solaiman & Anawar, 2015;Kushwaha et al., 2020). Therefore, it is of great significance to study the dynamics of rhizosphere microbiota in the cultivation of medicinal herbs.
The traditional Chinese medicinal herb "Beishashen" is the swelling root of perennial herb Glehnia littoralis Fr. Schmidt ex Miq (Shan & She, 1992). It contains a variety of coumarin compounds and alkaloids, and has a variety of activities such as antibacterial and anti-inflammatory (Taesook et al., 2010;Jing et al., 2021), so it is widely used in China and Southeast Asia. Due to the endangered status of the wild G. littoralis resources (Fu, 1992), the medicinal herb "Beishashen" on the market mainly comes from human cultivation. In China, G. littoralis has been cultivated for more than 600 years. Due to its special habitat requirements, it can only grow in sandy soil, which makes its planting area relatively fixed. Long-term planting in fixed areas results in typical negative plant-soil feedback (NPSF), leading to a decline in the yield and quality of the medicinal herbs (Cesarano et al., 2017;Yuan et al., 2022). Interactions between rhizosphere microbiota and plants affect plant health and productivity, and this process is very important in medicinal plant cultivation (Trivedi et al., 2020;Xiong et al., 2021). Given the importance of rhizosphere microbes, the study of rhizosphere microbiota may provide a way to improve the yield and quality of the medicinal herbs. Relevant microbial research has attracted the attention of the researchers, but so far, only the bacteriostatic activity of endophytic fungi and related studies have been reported (Huo et al., 2020).
To provide a basis for the interpretation and utilization of beneficial rhizosphere microbial resources, rhizosphere microbes of G. littoralis in genuine producing areas were analyzed by high-throughput sequencing technique. The composition, diversity, function, and dynamics of rhizosphere microorganisms at different development stages of G. littoralis, as well as the correlation between rhizosphere microorganisms and environmental factors, were investigated, hoping to provide reference data for optimizing cultivation of Chinese medicinal herb "Beishashen".

Sampling
All the samples were collected from Haiyang City (37 01′27.04″N, 120 44′53.71″E), Shandong Province, China, one of the main cultivation areas of G. littoralis. As a perennial plant, G. littoralis is usually harvested within a year when used as a Chinese medicinal herb. Therefore, samples growing for only one year in three growing seasons (Spring, Summer and Autumn) were collected according to the phenological stages, representing the seedling stage, the vigorous growth stage and the harvesting stage, respectively. In a 2 m × 2 m quadrat, three samples were collected on the diagonal of the quadrat as triplicates of each development stage. Since the root system at the seedling stage is small, in order to collect enough rhizosphere soil, the rhizosphere soil of 2-3 G. littoralis seedlings on the diagonal of each quadrat was collected as a replicate. The collected rhizosphere soil samples were quickly transported to the laboratory at low temperatures. Some soil samples were kept in the shade to test soil physiochemical properties, and the rest was stored at −80 C for high-throughput sequencing. Field experiments were approved by Ludong University (project number 20210301).

Measurement of soil physiochemical properties
Soil samples with roots and debris removed using a 2 mm sieve were air-dried and stored at 4 C for use. Fourteen common environmental factors in soil were detected according to previous studies, including available phosphorus (AP), available potassium (AK), ammonium nitrogen (AN), soil organic matter (SOM), total organic carbon (TOC) (Bao, 2000), nitrate nitrogen (NN) (Sun et al., 2016), Saccharase (SC), Urease (UE), alkaline phosphatase (AKP) (Guan, 1986); total nitrogen (TN), total hydrogen (TH), total carbon (TC), and total sulfur (TS) were detected by Elemental Analyzer (Elementar vario EL cube Elemental Analyzer; Elementar, Langenselbold, Germany) and the pH of soil samples was determined in 1:2.5 soil-water suspension using pH meter (Sartorius PB-10). The saccharase activity was expressed as mg glucose·d −1 ·g −1 soil, the urease activity was expressed as mg NH 3 -N·d −1 ·g −1 soil and the alkaline phosphatase activity was expressed as mg P 2 O 5 ·2h −1 ·g −1 soil.

DNA extraction and high-throughput sequencing
The total DNA of rhizosphere microorganisms was extracted using HiPure Soil DNA Kits (Tiangen, Beijing, China) according to the manufacturer's protocols. The extracted DNA was used for PCR amplification. For bacteria, the 16S V3-V4 region of the ribosomal RNA gene was amplified using primers 338F (5′-ACTCCTACGGGAGGCAGCA-3′) (Huse et al., 2008), and 806R (5′-GGACTACHVGGGTATCTAAT-3′) (Caporaso et al., 2011). The PCR program is as following: 95 C for 5 min, followed by 25 cycles at 95 C for 30 s, 50 C for 30 s, and 72 C for 40 s and a final extension at 72 C for 7 min. For fungi, the ITS (internal transcribed spacer) region of the ribosomal RNA gene was amplified using primers ITS-F (5′-CTTGGTCATTTAGAGGAAGTAA-3′) (Gardes & Bruns, 1993) and ITS-R (5′-GCTGCGTTCTTCATCGATGC-3′) (White et al., 1990) by PCR (95 C for 5 min, followed by 25 cycles at 95 C for 1 min, 50 C for 30 s, and 72 C for 1 min and a final extension at 72 C for 7 min). PCR products were confirmed by 1.8% agarose gel electrophoresis and purified using VAHTSTM DNA Clean Beads (Vazyme, Nanjing, China). A library was constructed, the Solexa PCR product was purified, and the sequenced library was constructed after quantification and homogenization. The DNA was purified with the Monarch DNA Gel Extraction Kit (Hongyue, Beijing, China) and then subjected to high-throughput sequencing on Illumina Novaseq 6000 (Biomarker Technologies, Beijing, China). The raw sequencing data has been uploaded to the NCBI Sequence Read Archive (SRA) database (BioProject: PRJNA903756). Three biological replicates of each stage were sequenced.

Feature of operational taxonomic units (OTUs)
The operational taxonomic units (OTUs) of rhizosphere bacteria and fungi were obtained from all the samples at the three growth stages (Fig. 1). A total of 1,885 bacterial OTUs were identified in rhizosphere soil of three growth stages, including 1,660 OTUs at the seedling stage, 1,267 OTUs at the vigorous stage, and 1,409 OTUs at the harvesting stage, accounting for 88.06%, 67.21% and 74.75% of the total bacterial OTUs at each growth stage, respectively. A total of 966 OTUs were shared by three different growth stages of G. littoralis, 244 OTUs specific to the seedling stage, 61 OTUs specific to the vigorous growth stage, and 95 OTUs specific to the harvesting stage ( Fig. 1 A). A total of 903 fungal OTUs were found in the soil at all growth stages, including 474 OTUs at the seedling stage

Variation in the abundance and diversity of rhizosphere bacteria and fungi
The bacterial OTUs were clustered into 25 phyla and 479 genera. At phylum level, the relative abundance of bacteria showed no obvious difference between the bacterial communities at the vigorous growth stage and harvesting stage, but obvious differences were found between the bacteria at the two stages and those at the seedling stage ( Fig. 2A). Proteobacteria was the major component of each bacterial community, representing 31.94%, 41.36% and 41.15% of the total species in each community at each development stage, respectively; Acidobacteriota was the second most abundant phylum after Proteobacteria, which happened in 22.30%, 19.61% and 17.94% of all the samples at each development stage, respectively; except for Proteobacteria and Acidobacteriota, Actinobacteria was the main phylum with abundance above 10%, which happened in 15.02%, 11.84% and 13.66% of the samples at each development stage, respectively (Table S1). The fungal OTUs were assigned to 10 phyla and 312 genera. Similar to bacteria community, the abundance of dominant phylum in fungi community showed no obvious difference between the vigorous growth stage and the harvesting stage, but there were obvious differences between the fungi at the two stages and those at the seedling stage ( Fig. 2A). Ascomycota was the most enriched phylum representing 70.64%, 78.98%, and 80.14% of the total species at each growth stage and Basidiomycota was the second abundant phylum after Ascomycota, which happened in 11.93%, 15.2% and 9.88% of the samples at each development stage (Table S2). Diversity of rhizosphere microbiota The alpha diversity of the bacterial community of G. littoralis showed significant differences between the seedling stage and the vigorous growth and harvesting stages, while there was no significant difference between the vigorous growth stage and the harvesting stage (Table 1, Fig. 3). The bacterial community's richness (ACE and Chao1) and diversity (Shannon index) were the highest at the seedling stage, and followed by the vigorous growth stage and the harvesting stage (Table 1). There was no significant difference in ACE, Chao1 and Shannon index between rhizosphere bacterial communities at the vigorous growth stage and the harvesting stage (Table 1). In addition, the richness and diversity of bacterial communities showed a gradually decreasing trend with the different development stages of G. littoralis (Table 1). The alpha diversity analysis results of fungi community are also listed in Table 1; it is shown that the fungi community' richness (ACE and Chao1) at the harvesting stage was the highest, followed by the seedling stage and the vigorous growth stage (Table 1). This result was different from that of the bacterial community. The Shannon index of fungi community at the seedling stage was the highest  and significantly higher than that at the vigorous growth stage. The Shannon index of fungi community at the harvesting stage was slightly higher than that at the vigorous growth stage, and the difference was small and not significant. PCoA analysis based on Bray-Curtis distances was performed to demonstrate changes in rhizosphere microbial community structure at different development stages. Both rhizosphere bacteria and rhizosphere fungi were divided into two groups (Fig. 4). One group included only microbiota at the seedling stages, and the other group included microbiota at the vigorous growth stage and the harvesting stage. This was consistent with the results of alpha diversity analysis, indicating that the rhizosphere microorganisms at the seedling stage are different from rhizosphere microorganisms at other development stages.

Characteristics of rhizosphere microbial function at different development stages
There was no significant difference in the abundance of the major function at class 2 KEGG pathways in the microbiota of three development stages (Fig. 5A). The dominant 10 functions of bacteria include "Global and overview maps", "Carbohydrate metabolism", "Amino acid metabolism", "Energy metabolism", "Metabolism of cofactors and vitamins", "Membrane transport", "Nucleotide metabolism", "Translation", "Signal transduction" and "Replication and repair" (Fig. 5A). "Global and overview maps" was the main component, accounting for 42.13-42.32% (Table S3).
Comparison analysis of the functions between the rhizosphere bacteria of G. littoralis at three development stages showed that six, three and seven levels of functions were significantly different between the bacteria at the seedling stage and the vigorous growth stage, the vigorous growth stage and the harvesting stage, and the seedling stage and the harvesting stage, respectively (Figs. 5C-5E). From the seedling stage to the vigorous growth stage, the function of "Circulatory system" and "Endocrine and metabolic diseases" increased significantly, while "Glycan biosynthesis and metabolism", "Transport and catabolism", "Nervous system" and "Substance dependence" showed a significant decrease (Fig. 5C). Comparison of functions between the vigorous growth stage and the harvesting stage (Fig. 5D) showed a marked increase in "Replication and repair", while "Drug resistance" and "Cancers" presented a marked decline. From the vigorous growth stage to the harvesting stage (Fig. 5E), the functions of "Cell growth and death", "Carbohydrate metabolism", "Infectious diseases: parastitic" and "Substance dependence" significantly decreased, and the functions of "Circulatory system", "Immune system" and "Environmental adaptation" significantly increased. Eight trophic modes were found in the fungi at three development stages (Fig. 5B), including symbiotroph, saprotroph, pathotroph, saprotroph-symbiotroph, pathotroph-symbiotroph, pathotroph-saprotroph, pathotroph-saprotroph-symbiotroph, and pathogen-saprotroph-symbiotroph. Saprotroph was the main trophic mode, accounting for 39.30-60.59% of the total fungi at each stage, and its relative abundance gradually increased with the development stages of G. littoralis (Fig. 5B, Table S4). In addition, pathotroph-saprotroph-symbiotroph was the second trophic mode after Saprotroph, accounting for 7.11-16.36% of fungi at each stage (Table S4).
Functional differences analysis showed that only the abundance of pathotrophsaprotroph-symbiotroph were significantly different at different development stages, and revealed an increasing trend and then a decreasing trend with the development of G. littoralis (Figs. 5F and 5G). Finally, the PCA analysis of the functions at the three development stages was performed, and distinct differentiation was observed between the functions at the seedling stage and those at the vigorous growth and harvesting stages (Fig. S1).

Microbial biomarkers in the rhizosphere at different development stages
The random forest method was used to identify microbial biomarkers in the rhizosphere of G. littoralis. The optimal model of rhizosphere bacteria and fungi (Figs. 6A and 6D) was established based on the cross-validation results. The biomarkers identified in the samples were sorted according to their importance to the community (Mean Decrease Accuracy, MDA) from large to small to obtain a histogram (Figs. 6B and 6E). A heatmap analysis of rhizosphere microbes at different development stages was also performed (Figs. 6C and These biomarkers can be used to distinguish samples at different stages accurately (Fig. 6C). The five orders with the greatest impact on bacterial community were Bdellovibrionales, Actinomycetales, Erysipelotrichales, Enterobacterales, and Lactobacillales (Fig. 6B). At fungal order level, 22 important biomarkers were identified (Figs. 6D and 6E). These biomarkers were also able to accurately distinguish samples at different development stages. The top five orders were Mucorales, Cystobasidiales, Ustilaginales, GS11, Wallemiales (Fig. 6E). The clustering results showed that both the bacterial and fungal biomarkers at three development stages clustered into two branches (Figs. 6C and 6F), one was composed with the microbiota at the seedling stage, another was composed with the microbiota at the vigorous growth stage and at the harvesting stage, which was consistent with the results of the above PcoA analysis (Figs. 4A and 4B).

The interactions of microbes-microbes and microbes-environmental factors at different development stages
Redundancy analysis (RDA) revealed that eight environmental factors drove changes in the microbiome at different development stages (Figs. 7A and 7C). In the rhizosphere bacterial community, soil organic matter (SOM), pH, urease (UE), and available phosphorus (AP) showed a positive correlation with the bacterial community at the seedling stage; saccharase (SC) and nitrate nitrogen (NN) were positively correlated with the bacterial community at the vigorous growth stage and the harvesting stage; alkaline phosphatase (AKP) and available potassium (AK) impacted the bacterial community at all stages (Fig. 7A).
The bacteria that played a major role at the seedling stage were Acidobacteriota, Actinobacteriota, Armatimonadota, Bacteroidota, Chloroflexi, Dependentiae, Firmicutes, Fusobacteriota, Gemmatimonadota, Methylomirabilota, Nitrospirota, Planctomycetota, and Verrucomicrobiota; The main bacteria at the vigorous growth and the harvesting stages were Proteobacteria, Patescibacteria, Bdellovibrionota, Campylobacterota, Cyanobacteria, Desulfobacterota, Elusimicrobiota, and Myxococcota (Fig. 7B). In the rhizosphere fungal community, SOM, AP, pH, AK, AKP and UE were the main influencing factors at the seedling stage, while SC and NN were the main influencing factors at the vigorous growth and harvesting stages (Fig. 7C). Chytridiomycota, Mortierellomycota, Olpidiomycota were the main compositions of fungal community at the seedling stage, and Ascomycota, Basidiomycota, Glomeromycota, Mucoromycota, and Rozellomycota were the main compositions of fungal community at the vigorous growth and harvesting stages (Fig. 7D). Network analysis was carried out to understand interactions between different rhizosphere microbiota, and between rhizosphere microbiota and environmental factors. It was shown that rhizosphere bacterial community had more complex interaction network than rhizosphere fungal community (Figs. 8A and 8B), which may be attributed to the higher abundance and diversity of the bacterial community. In bacterial community, all 13 environmental factors revealed positive or negative impacts on bacteria at the Phylum level, with more negative (56.60%) than positive correlations (43.40%) were observed (Fig. 8A). Planctomycetota contributed the most to the network, followed by Firmicutes, Myxococcota, and Proteobacteria. In fungal community, only 7 of the 13 environmental factors were associated with rhizosphere fungi (Fig. 8B). AN may be the key factors driving the composition of rhizosphere fungal community. Chytridiomycota, Mortierellomycota, and Olpidiomycota may be the three phyla that had the greatest influences on the fungal community (Fig. 8B). Interestingly, total hydrogen (TH) and total sulfur (TS) were negatively correlated with rhizosphere bacteria but positively correlated with rhizosphere fungi, suggesting that rhizosphere bacteria and fungi respond differently to environmental factors. A co-occurrence network of bacteria-fungi-environmental factors was established (Fig. S2), and more positive interactions (60%) were found in bacteria-bacteria and bacteria-fungi. The number of OTUs in both bacteria and fungi communities first decreased and then increased along the seedling-vigorous growth-harvesting development stages. In addition, the proportion of OTUs special to each stage of bacterial and fungal communities also showed similar results, with the values of 14.70-4.81-6.74% and 33.33-20.68-41.46%, respectively. The results may reflect the process of rhizosphere microorganisms colonization. At the seedling stage, due to the weak root system of G. littoralis, OTUs may represent indigenous microorganisms in the soil. At the harvesting stage, the rhizosphere microbes adapted to the root exudates are gradually enriched, and community structure become more stable, thus total and special OTUs number increases. In terms of species abundance, Proteobacteria, Acidobacteria, and Actinobacteria are the main components of bacterial community, and Ascomycota and Basidiomycota are the dominant fungi at each development stage. This is similar to the species composition of the rhizosphere microorganisms in other plants (Nie et al., 2012;Li et al., 2020). The composition of rhizosphere bacteria of G. littoralis is also consistent with the endophytic bacterial composition previously reported (Huo et al., 2020). This is as expected, as previous studies have confirmed that most endophytic bacteria come from the soil (Huo et al., 2020). There have been no reports on the composition of endophytic fungi of G. littoralis. Proteobacteria play an important role in plant growth and development by promoting photosynthesis to improve the utilization of C-source (Ma et al., 2021). The abundance of Proteobacteria increased significantly from the seedling stage to the vigorous growth stage, but remained stable from the vigorous growth stage to the harvesting stage, which may be related to changes in soil C and N content. Increased abundance of Proteobacteria leads to more efficient use of complex carbohydrates and promotes plant adaptation to high C and N environments. The abundance of Acidobacteria decreased gradually from the seedling stage to the harvesting stage. Acidobacteria mainly degraded plant residues, participated in the metabolism of carbon compounds, and photosynthesis. The reason for the decline of abundance in the late development of G. littoralis may be that the microbial community used simple amino acids in the early stage and complex carbohydrates in the late stage. These characteristics have a positive impact on the growth of G. littoralis. Ascomycota is the largest fungal group and is important to rhizosphere ecosystems (Beimforde et al., 2014). This group contains a large number of beneficial species (Braun et al., 2012) and harmful species (Peterson et al., 2021). Basidiomycota is an important pathogen lineage in the fungal kingdom (Coelho et al., 2017), and is the main fungi responsible for rust disease of G. littoralis. It was relatively abundant at three development stages (9.88-15.2%), with the greatest abundance occurring at the vigorous growth stage, suggesting that the invasion of pathogenic fungi may have occurred at the vigorous growth stage, but G. littoralis remained healthy. Previous studies have shown that Proteobacteria and Actinobacteria are important phylum associated with plant disease suppression (Mendes et al., 2011). Therefore, it is speculated that Proteobacteria and Actinobacteria with high abundance may participate in the disease suppression of G. littoralis and thus maintain the health of G. littoralis. Compared with the previously reported microbial abundance in soils at the harvesting stage with that in continuous crop obstacles soil of G. littoralis (Gu et al., 2021), the abundance of Proteobacteria (41.15% vs 33%) and Basidiomycota (9.88% vs 3.8%) was significantly higher, while that of Actinobacteria (13.66% vs 13%) and Acidobacteria (17.94% vs 17%) was similar to the former. These results suggest that Proteobacteria with higher abundance may play a key role in disease suppression. In cultivation of G. littoralis, controlling the occurrence of diseases and pests is an important prerequisite to ensure the yield. Grasping the occurrence time and abundance of main phytobacteria can provide reference for field management and yield increase at development stages of G. littoralis. The alpha diversity of bacterial community showed a decline with the seedling-vigorous growth-harvesting stage, consistent with results in other species (Shi et al., 2016;Cordovez et al., 2019). The alpha diversity of fungal community decreased first and then increased with the development time. These results indicate that soil type is changing from a bacteria-dominated state to a fungi-dominated stage. However, the Ace, Chao1 and Shannon index of bacteria were significantly higher than those of fungi at all development stages. This suggests that the bacterial species in the rhizosphere soil are more abundant and diverse than fungi. Previous studies have suggested that bacteria-dominated soils are healthier. As the composition of rhizosphere soil changes from a bacteria-dominated state to a fungi-dominated stage, diseases and pests increase, leading to lower crop yield and other adverse effects (Liao et al., 2011). The results presented here indicate that the healthy state of the rhizosphere may be the basis to ensure the high quality of the medicinal herbs of "Beishashen" at the harvesting stage. Combined with the changes in microbial abundance, the difference in alpha diversity index between bacteria and fungi may represent a change trend in soil health status in the late development stage of G. littoralis, and the high abundance of Proteobacteria maintained the healthy development of G. littoralis.
High functional diversity was observed in the rhizosphere bacterial and fungal communities of G. littoralis. For bacterial community, the 10 key functions dominate at all the development stages. Of these, the predictive function of "Membrane transport" was thought to be associated with symbiotic interactions between bacterial community and other organisms (Burke et al., 2011;Ma et al., 2021), suggesting that rhizosphere bacteria may mediate rhizosphere interactions throughout the whole growth process of G. littoralis. The functions of "Carbohydrate metabolism", "Energy metabolism", and "Amino acid metabolism" are associated with carbon and nitrogen fixation, and rhizosphere bacteria are speculated to convert carbon and nitrogen from soil, thus providing available substances for plants. Through the functional difference analysis, significant differentiations of the functions at the early development stage and the middle to late development stages were observed, which is supported by the results of PCA analysis. The abundance of "environmental adaptation" function at the harvesting stage was significantly higher than that at the seedling stage. This function mainly related to "plant-pathogen interaction", indicating that the bacterial community may be involved in plant disease resistance at the harvesting stage. For fungal community, functional differentiation also occurred between the early development stage and the middle to late development stages. The relative abundance of saprotroph, symbiotroph, and pathotroph-saprotroph-symbiotroph fungi increased gradually with the growth and development of G. littoralis. These fungi may accelerate the decomposition of soil organic matter and convert insoluble minerals in the soil into nutrients available to G. littoralis. The results of functional difference analysis also confirmed functional differentiation between early and late development: Only functional changes were found at the seedling stage and vigorous growth stage, and at the vigorous growth stage and harvesting stage, but not at the seedling stage and harvesting stage.
Biomarkers at three development stages were identified using a random forest machine learning algorithm. The 47 bacterial and 22 fungal biomarkers exhibit different correlations at different developmental stages, and these characteristics can be used to distinguish G. littoralis at different development stages. Biomarkers may be useful indicators for identifing plant origin, and some studies propose using biomarkers to determine the origin of imported soybeans (Zhou et al., 2022). As a traditional Chinese medicine, the quality of "Beishashen" is of paramount importance. Medicinal herbs grown in a specific area have better efficacy and are called "authentic medicinal herbs". The samples collected in this study are from the authentic producing areas of G. littoralis, therefore this study can provide an option for reliable medicinal herbs origin traceability technology.
The structure and diversity of rhizosphere microorganisms are influenced by multiple environmental factors. The RDA analysis showed that eight environmental factors had significantly influence on the rhizosphere microbial community, and different environmental factors had different influence on microbial community at different development stages. SOM, pH, UE, and AP positively impacted the bacterial and fungal communities at seedling stage, among which the influence of pH on bacterial and fungal communities was the most remarkable, and bacterial community was significantly affected by SOM, while fungal community was significantly impacted by AP. These results demonstrate the importance of pH, SOM and AP in early development of G. littoralis. Previous studies have shown that rhizosphere microbial communities promote colonization of beneficial microbes by regulating pH to suppress plant immunity (Yu et al., 2019). This suggests that pH may play an important role in rhizosphere microorganisms colonization at the seedling stage of G. littoralis. SOM and AP are important indicators of soil fertility, and their effects on early plant rhizosphere need further study. SC and NN showed significant positive effects on bacterial and fungal communities at the vigorous growth stage and harvesting stage. SC can hydrolyse sucrose in soil to produce small molecules of glucose, which is an important carbon source for most microbes. NN is an important nitrogen source for most microorganisms. Therefore, it can be inferred that rhizosphere microbes absorbs a lot of carbon and nitrogen sources and carried out vigorous metabolic activities at the vigorous growth stage and harvesting stage of G. littoralis. The metabolites produced can provide available nutrients for G. littoralis and promote the production of active secondary metabolites. RDA analysis also provides a list of the bacteria and fungi that play a positive role in different development stages, so the microorganisms that are beneficial to the growth of G. littoralis can be screened. Isolation and culture of these beneficial microorganisms are of great significance to improve the yield and quality of G. littoralis. Understanding the influence of various environmental factors on rhizosphere microbial community can help us to construct beneficial microbiome by regulating environmental factors to promote the healthy growth of plants and improve the yield of medicinal herbs.
Co-occurrence networks can identify putative interactions between microorganisms in the environment (Berry & Widder, 2014). The results revealed more complex co-occurrence networks in bacterial community than in fungal community. The formation of bacterial community involves more environmental factors than the formation of fungal community. Other studies have found complex co-occurrence networks in rhizosphere bacterial communities (Fan et al., 2017;Zhang et al., 2022). Some studies have shown that rhizosphere influences ultimately lead to a decrease in bacterial community diversity and more complex symbiotic networks (Shi et al., 2016;Cordovez et al., 2019). Our findings support this conclusion. In order to understand all rhizosphere microbial interactions, a co-occurrence network of bacteria-fungi-environmental factors was established. More positive interactions (60%) in bacteria-bacteria and bacteria-fungi were found. The results suggest that the rhizosphere microorganisms form a complex mutualistic symbiosis network, which might be beneficial to the growth and development of G. littoralis. These results are also consistent with previous studies showing that most bacteria cluster together as functional groups, which use plant-derived resources more efficiently and provide a greater number of services (Saleem, Hu & Jousset, 2019).

CONCLUSIONS
In this study, high-throughput sequencing technology revealed the dynamic changes of rhizosphere microbial communities (bacteria and fungi) at different development stages of G. littoralis. The composition, diversity and function of rhizosphere bacterial and fungal communities are closely related to the development stage of G. littoralis. Eight environmental factors play a vital role in driving rhizosphere microbial changes at different development stages. This study provides data support for understanding the structure and composition of the rhizosphere microbial community in the development period of G. littoralis, and also lays a foundation for improving the yield and quality of G. littoralis by regulating the microbial community in the future. There are still many questions to be explored about how to utilize rhizosphere microbial resources, how to increase microbial diversity in agro-ecosystems, and the mechanisms of microbial diversity contribute to agricultural production.

ADDITIONAL INFORMATION AND DECLARATIONS Funding
This work was supported by the National Natural Science Foundation of China (Grant Number 32070355), the Shandong Provincial Natural Science Foundation (Grant Number ZR2020MC030), and the Shandong Provincial Agricultural Elite Varieties Project (Grant Number 2019LZGC01805). There was no additional external funding received for this study. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Ailan Wang conceived and designed the experiments, analyzed the data, prepared figures and/or tables, authored or reviewed drafts of the article, and approved the final draft.

Field Study Permissions
The following information was supplied relating to field study approvals (i.e., approving body and any reference numbers): Field experiments were approved by the Ludong University (project number 20210301).